#include "ZigguratRandomRoutine.h"

static unsigned long jz,jsr=123456789;

#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define UNI (.5 + (signed) SHR3*.2328306e-9)
#define IUNI SHR3

static long hz;
static unsigned long iz, kn[128], ke[256];
static float wn[128],fn[128], we[256],fe[256];

#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())

/* nfix() generates variates from the residue when rejection in RNOR occurs. */

float SmartMathLibrary::InternalRoutines::ZigguratRandomRoutine::Nfix(void)
{
	const float r = 3.442620f;     /* The start of the right tail */
	static float x, y;
	for(;;)
	{  x=hz*wn[iz];      /* iz==0, handles the base strip */
	if(iz==0)
	{ do{ x=-System::Math::Log(UNI)*0.2904764; y=-System::Math::Log(UNI);}	/* .2904764 is 1/r */
	while(y+y<x*x);
	return (hz>0)? r+x : -r-x;
	}
	/* iz>0, handle the wedges of other strips */
	if( fn[iz]+UNI*(fn[iz-1]-fn[iz]) < System::Math::Exp(-.5*x*x) ) return x;

	/* initiate, try to exit for(;;) for loop*/
	hz=SHR3;
	iz=hz&127;
	if(System::Math::Abs(hz)<kn[iz]) return (hz*wn[iz]);

	}

}

/* efix() generates variates from the residue when rejection in REXP occurs. */
float SmartMathLibrary::InternalRoutines::ZigguratRandomRoutine::Efix(void)
{ float x;
for(;;)
{  if(iz==0) return (7.69711-System::Math::Log(UNI));          /* iz==0 */
x=jz*we[iz]; if( fe[iz]+UNI*(fe[iz-1]-fe[iz]) < System::Math::Exp(-x) ) return (x);

/* initiate, try to exit for(;;) loop */
jz=SHR3;
iz=(jz&255);
if(jz<ke[iz]) return (jz*we[iz]);
}
}
/*--------This procedure sets the seed and creates the tables------*/

void SmartMathLibrary::InternalRoutines::ZigguratRandomRoutine::Zigset(unsigned long jsrseed)
{ 
const double m1 = 2147483648.0, m2 = 4294967296.;
double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
int i;
jsr^=jsrseed;

/* Set up tables for RNOR */
q=vn/System::Math::Exp(-.5*dn*dn);
kn[0]=(dn/q)*m1;
kn[1]=0;

wn[0]=q/m1;
wn[127]=dn/m1;

fn[0]=1.;
fn[127]=System::Math::Exp(-.5*dn*dn);

for(i=126;i>=1;i--)
{dn=System::Math::Sqrt(-2.*System::Math::Log(vn/dn+System::Math::Exp(-.5*dn*dn)));
kn[i+1]=(dn/tn)*m1;
tn=dn;
fn[i]=System::Math::Exp(-.5*dn*dn);
wn[i]=dn/m1;
}

/* Set up tables for REXP */
q = ve/System::Math::Exp(-de);
ke[0]=(de/q)*m2;
ke[1]=0;

we[0]=q/m2;
we[255]=de/m2;

fe[0]=1.;
fe[255]=System::Math::Exp(-de);

for(i=254;i>=1;i--)
{de=-System::Math::Log(ve/de+System::Math::Exp(-de));
ke[i+1]= (de/te)*m2;
te=de;
fe[i]=System::Math::Exp(-de);
we[i]=de/m2;
}
}